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CNj , This paper gives an introduction to the Keldysh formalism, with emphasis on its usefulness in 

time- dependent density functional theory. In the first part we introduce the Keldysh contour and the 
one-particle Green function defined on this contour. We then discuss how to combine and manipulate 
functions with time-arguments on the contour. The effects of electron-electron interaction can be 
taken systematically into account, as we illustrate by propagating the Kadanoff-Baym equations for 
the second-order self-energy approximation. It is important to use conserving approximations such 
that the evolution of the electron density, momentum, and total energy agrees with the macroscopic 
conservation laws. One of the main topics in this paper is the non-interacting Green function, which 
is the relevant quantity for time-dependent density functional theory. We discuss this Green function 
in detail, and show how the Keldysh contour in a simple way allows us to derive the time-dependent 
Kohn-Sham potential from an action functional. The formalism in a similar way leads to response 
functions that obey the causality principle. To illustrate these points, we discuss the time-dependent 
optimized effective potential equations. 
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I. INTRODUCTION 



We will in this paper give an introduction to the Keldysh formalism, which is an extremely useful tool for first- 
principles studies of nonequilibrium many-particle systems. Of particular interest for TDDFT is the relation to 
non-equilibrium Green functions (NEGF), which allows us to construct exchange-correlation potentials with memory 
by using diagrammatic techniques. For many problems, such as, e.g., quantum transport or atoms in intense laser 
pulses, one needs exchange-correlation functionals with memory, and Green function techniques offer a systematic 
method for developing these. The Keldysh formalism is also necessary for defining response functions in TDDFT and 
for defining an action functional needed for deriving TDDFT from a variational principle. We will in this section 
give an introduction to the nonequilibrium Green function formalism, intended to illustrate the usefulness of the 
theory The formalism does not differ much from ordinary equilibrium theory, the main difference being that all 
, time-dependent functions are definied for time-arguments on a contour, known as the Keldysh contour. 

The Green function, G(f, t; r 1 , t') is a function of two space- and time-coordinates, and is obviously more complicated 
I/"") | than the one-particle density n(f,t), which is the main ingredient in TDDFT. However, the advantage of the NEGF 
methods is that we can systematically improve the approximations by taking into account particular physical processes 
(represented in the form of Feynman diagrams) that we believe to be important. The Green function provides us 
directly with all expectation values of one-body operators (such as the density and the current), and also the total 
energy, ionization potentials, response functions, spectral functions, etc.. In relation to TDDFT, this is useful not only 
for developing orbital functionals and exchange-correlation functionals with memory, but also for providing insight in 
(— i ■ the exact properties of the non-interacting Kohn-Sham system. 

In the following, we shall focus on systems that are initially in thermal equilibrium. We will start by introducing the 
Keldysh contour and the nonequilbrium Green functions, and then explain how to combine and manipulate functions 
with time variables on the contour. While we in TDDFT take exchange- and correlation-effects into account through 
•i-H , u xc [n], the corresponding quantity in Green function theory is the self-energy E[G]. Just like v xc , the self-energy 
functional must be approximated. For a given functional E[G], it is important that the resulting observables obey the 
macroscopic conservation laws, such as, e.g., the continuity equation. These approximations are known as conserving, 
and will be discussed briefly In the last part of this section we will discuss the applications of the Keldysh formalism 
in TDDFT, including the relation between E and v xc , the derivation of the Kohn-Sham equations from an action 
functional, and the derivation of an f KC functional. As an illustrative example, we will discuss the time-dependent 
exchange-only optimized effective potential approximation. 
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II. THE KELDYSH CONTOUR 



In quantum mechanics we associate with any observable quantity O a hermitean operator O. The expectation value 
Tr {pqO} gives the value of O when the system is described by the density operator p and the trace denotes a sum 
over a complete set of states in Hilbert space. For an isolated system the Hamiltonian Hq does not depend on time, 
and the expectation value of any observable quantity is constant, provided [pq,Hq] = 0. In these notes we want to 
discuss how to describe systems that are isolated for times t < 0, such that H(t < 0) = Hq, but disturbed by an 
external time-dependent field at t > 0. The expectation value of O at t > is then given by the average on the initial 
density operator po of the operator O in the Heisenberg representation, 



0{t) = (d H (t))=Ti-{pod H (t)} = TY{poS(0;t)dS(t;0)}, 

where the operator in the Heisenberg picture has a time-dependence according to 0/f(<) 
evolution operator S(t;t') is the solution of the equations 

»4-5(f; t') = H(t)S(t; t') and i-^-;S(t; t') = -S(t; t')H(t'), 
at at' 

with the boundary condition S(t;t) = 1. It can be formally written as 



S(t;t') = 



Te -if*diH(.i) t>t , 
Te -ij t UiH(F) t<t , 



(1) 

§(0;t)6S(t;0). The 
(2) 

(3) 



In Eq. ©, T is the time-ordering operator that rearranges the operators in chronological order with later times 
to the left; T is the anti-chronological time-ordering operator. The evolution operator satisfies the group property 
S(t;ti)S(ti;t') — S(t;t') for any t\. Notice that if the Hamiltonian is time-independent in the interval between t 

and t' , then the evolution operator becomes S(t;t') = e^ 111 ^^ 1 \ If we now let the system be initially in thermal 
equilibrium, with an inverse temperature f3 = 1/ksT and chemical potential p,, the initial density matrix is po = 
g-^(J?o-/iW)y"ji r | e -^(J?o-/iJV)j._ Assuming that Ho and iV commute, po can be rewritten using the evolution operator 

S with a complex time- argument, t — —i(3, according to p Q = e^^ N S{-i(3; 0)/Tr {e^ A '5(-i/3; 0)}. Inserting this 
expression in (IJ, we find 



0(t) = 



Tr 


epi**§(-i0 ; 0)5(0; t)OS(t;0) 


Tr 


ePrf§(-i0\ 0) 





(4) 



Reading the arguments in the numerator from the right to the left, we see that we can design a time-contour 7 with 
a forward branch going from to t, a backward branch coming back from t and ending in 0, and a branch along the 
imaginary time-axis from to — ij3. This contour is illustrated in Fig. ^ Note that the group property of S means 
that we are free to extend this contour up to infinity. We can now generalize Q , and let z be a time-contour variable 
on 7. Letting the variable z run along this same contour, (@J) can be formally recast as 



O(z) 



Tr 




Tr 







(5) 



The contour ordering operator T c moves the operators with "later" contour variable to the left. In JSJ), 0{z) is not 
the operator in the Heisenberg representation [the latter is denoted with Off(t)]. The contour-time argument in O 
is there only to specify the position of the operator O on 7. A point on the real axis can be either on the forward 
(we denote these points or on the backward branch (denoted t + ), and a point which is earlier in real time, can 
therefore be later on the contour, as illustrated in Fig. 

If z lies on the vertical track, then there is no need to extend the contour along the real axis. Instead, we have 



0{z) 



Tr 




f„-"d**od e - 


-if*dz6 ~ 


Tr 


e -0(H o -nN) 





Tr 



e -/3(H -MJV)O 



Tr e -P(Ho-nN) 
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FIG. 1: The Keldysh contour, starting at t = 0, and ending at t = — i/3, with f on the backward branch and t' on the forward 
branch. By definition, any point lying on the vertical track is later than a point lying on the forward or backward branch. 

where the cyclic property of the trace has been used. The right hand side is independent of z and coincides with the 
thermal average Tr{poO}. It is easy to verify that JSJ) would give exactly the same result for 0(t), where t is real, if 
the Hamiltonian was time- independent, i.e. H (t) = Hq also for t > 0. 

To summarize, in the variable z lies on the contour of Fig. ^ the- r.h.s. gives the time-dependent statistical 
average of the observable O when z lies on the forward or backward branch, and the statistical average before the 
system is disturbed when z lies on the vertical track. 

III. NONEQUILIBRIUM GREEN FUNCTIONS 

We now introduce the nonequilibrium Green function (NEGF), which is a function of two contour time- variables. 
In order to keep the notation as light as possible, we here discard the spin degree of freedom; the spin index may 
be restored later as needed. The field operators ip(r), ^y{f) destroy and create an electron in r and obey the 
anticommutation relations {ip(r), ^' (r*)} = d(f— ?). We write the Hamiltonian H(t) as the sum of a quadratic term 



h(t) = / dfdr^ ^(r)(r\h(t)\r')ip(r') 



and the interaction operator 



(6) 



(7) 



We use boldface to indicate matrices in one-electron labels, e.g., h is a matrix and (r^hlr*) is the (r, f) matrix 
element of h. When describing electrons in an electro-magnetic field, the quadratic term is given by (r|h(i)|r') = 

5{r - f 1 ) \[V/i + A{r, t)} 2 /2 + i;(r,t)}. 

The definition of an expectation value in Q can be generalized to the expectation value of two operators. The 
Green function is defined as 



(8) 



where the contour variable in the field operators specifies the position in the contour ordering. The operators have 
a time-dependence according to the definition of the Heisenberg picture, e.g. ^^(r, z) = S(0; z)^(r)S(z;0). Notice 
that if the time-argument z is located on the real axis, then i/)jy(r, < + ) = tpn {r, t-). If the time-argument is on the 
imaginary axis, then ip(r, —it) is not the adjoint of 4>(r, — ir) since (~ir; 0) ^ S'(0; —it). The Green function can 
be written 



GO; z') = 0(z, z')G > (z; z') + 6{z', z)G < (z; z'). 



(9) 
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The function 6(z, z') is denned to be 1 if z is later on the contour than z', and otherwise. From the definition of the 
time-dependent expectation value in Eq. Q, it follows that the greater Green function G > (z; z'), where z is later on 
the contour than z' is 



G > (r,z;r,z') = 

If z' is later on the contour than z, then the Green function equals 

G < (r,z;r,z') = -- 



Tr 






Tr 


ePrf§(-%P; 0) 





(10) 



Tr 




-if3;Q)^ H {?,z')i> H {r,z) 


Tr 


e PnNS(-i/3;0) 





(11) 



The extra minus sign on the right hand side comes from the contour ordering. More generally, rearranging the 
field operators and ifi (later arguments to the left), we also have to multiply by (— l) p , where P is the parity 
of the permutation. From the definition of the Green function, it is easily seen that the electron density, n(r, z) — 

{ty[i(r, z)%pH (r, z)) and current is obtained according to 



n(r, z) = —iG(r, z; r, z + ) 



(12) 



j(f,z) = -i 



V 



2i 2i 



- + A(r, z) 



G(?,z;7,z') 



(13) 



where z + indicates that this time-argument is infinitcsimally later on the contour. 

The Green function G(z; z') obeys an important cyclic relation on the Keldysh contour. Choosing z = 0_, which 
is the earliest time on the contour, we find G(0_;z') = G < (0;z'), given by Hll|) with t/jjiir, 0) = ip(r). Inside the 
trace we can move ip(r) to the left. Furthermore, we can exchange the position of ?/>(r) and e^^ N by noting that 
■4>(r)e^ = e^+^tp(f). Using the group identity §(-ip; 0)5(0; -iP) = 1, we obtain 



G(r,0-;7,z') 



l Tr 


'i> H (7)e^S(-ip; 0)^(^,2/)' 


i 


Tr 


eP^S(-ip;0) 





e ^Tr 


eP»*S{-ip\ 0)^ H (r, -ip)ft H (f, z') 




e HnNS{-ip] 0) 





(14) 



The r.h.s. equals —e^ tJl (r\G(—iP;z')\f / ). Together with a similar analysis for G(z;0_), we conclude that 

G(0_;z') = -e^G(-ip;z') and G(z; 0_) = -e- p »G{z\ -i0). (15) 

These equations constitute the so called Kubo-Martin-Schwinger (KMS) boundary conditions 0,0- From the defini- 
tion of the Green function in iJHJ, it is easily seen that the G(z; z) has a discontinuity in z — z', 



G>{z;z) = G<{z-z)-i\. 



(16) 



Furthermore, for both time-arguments on the real axis we have the important symmetry [G^(£';t)] = — G^(t;t'). 
As we shall see, these relations play a crucial role in solving the equation of motion. 

IV. THE KELDYSH BOOK-KEEPING 

The Green function belongs to a larger class of functions of two time-contour variables that we will refer to as 
Keldysh space. These functions can be written on the form 



k(z; z') = S(z, z')k d (z) + 0(z, z')k>(z- z') + 6(z', z)k<(z; z'), 



(17) 
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where the <5-function on the contour is defined as S(z, z') — d6(z, z')/dz 21]. These functions are somewhat compli- 
cated due to the fact that each of the time-arguments can be located on three different brances of the contour, as 
illustrated in Fig. ^ Below we systematically derive a set of identities that are commonly used for dealing with such 
functions and will be used extensively in the following sections. Most of the relations are well known while others, 
equally important 5|, are not. Our aim is to provide a self-contained derivation of all of them. A table at the end 
of the Section summarizes the main results. For those who are not familiar with the Keldysh contour, we strongly 
recommend to scan what follows with pencil and paper. 

It is straightforward to show that if a(z; z') and b(z; z') belong to the Keldysh space, then 

c(z;z') = I dz a(z; z)b{z; z') (18) 

also belongs to the Keldysh space. For any k(z; z') in the Keldysh space we define the greater and lesser functions on 
the physical time axis 

k>(t;t') = fc(t+;t'_), k<(t,t') = fc(t_;i+). 

We also define the following two-point functions with one argument t on the physical time axis and the other r on 
the vertical track 

k (t; r) ee k(t±; r), fe r (r, t) = fc(r; t±). (19) 

In the definition of fcl and we can arbitrarily choose t + or t_ since r is later than both of them. The symbols 
"]" and "["" have been chosen in order to help the visualization of the time arguments. For instance, "]" has a 
horizontal segment followed by a vertical one; correspondingly, fcl has a first argument which is real (and thus lies on 
the horizontal axis) and a second argument which is imaginary (and thus lies on the vertical axis). We will also use 
the convention of denoting the real time with latin letters and the imaginary time with greek letters. 

If we write out the contour integral in l|18|) in detail, we see with the help of Fig. ^ that the integral consists of 
four main parts. First, we must integrate along the real axis from z — 0_ to z = t'_, for which a = a > and b = b < . 
Then, the integral goes from z = t'_ to z ~ f+, where a — a> and b = b y . The third part of the integral goes along 
the real axis from z — t+ to z — 0+, with a = a < and b = b > . The last integral is along the imaginary track, from 
0+ to — i(3, where a = a} and b = b^ . In addition, we have the contribution from the singular parts, a s and b s , which 
is trivial since these integrals involve a <5-function. With these specifications, we can drop the i-subscripts on the 
time-arguments and write 

c> (*;*') = a > (t,t')b s (t')+a s (t)b > {t,t')+ f dt a> (t;t)b<(i;t') 

Jo 

+ / dta > (t;t)b > (t;t')+ / dt a < (t;t)b > (t;f) + / d? a) (t;f)tf (?;t'). 
Jt' Jt Jo 

The second integral on the r.h.s. is an ordinary integral on the real axis of two well defined functions and may be 
rewritten as 

t rO r-t 

dt a> (t; t)b> (t; t') = / dt a> (t; t)b> (t; t') + / dt a> (t; t)b> (i; t'). 
It' Jt' Jo 

Using this relation, the expression for c 5 " becomes 

c>(t;t') = a > (t,t')b s (t')+a s (t)b > (t,t')- f dt a> (t; t) [b> (f; t') - b< (i; t')] 

Jo 

ft r-i/3 

+ / dt[a > (t;t)-a < (t;t)]b > (t;t')+ df a 1 (t; f )6 r (r; t'). (20) 



jo Jo 
Next, we introduce two other functions on the physical time axis 



k R (t;t') ee S(t,t')k 5 + e(t~t')[k > (t;t') -k<{t-t')], (21) 
k A (t;t') ee 5(t,t')k s -6(t' -t^frt^-^frt')}. (22) 
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The retarded function k R (t;t') vanishes for t < t', while the advanced function k A (t;t') vanishes for t > t'. The 
retarded and advanced functions can be used to rewrite (|2U[I in a more compact form 

poo r — i0 

c>(t;t')= / di[a > (t;t)b A (i;t') + a R (t;i)b > (t-t r )}+ / df a 1 (t; f)6 r (f; t'). 



It is convenient to introduce a short hand notation for integrals along the physical time axis and for those between 
and —i(3. The symbol "•" will be used to write J °° dtf(t)g(t) as / • g, while the symbol "★" will be used to write 

df/(f) ff (f) as/* 5 . Then 



c> =a> -b A +a R -b> +aUtf. (23) 
Similarly, one can prove that 

c< =a< -b A +a R -b< +aUtf. (24) 
Equations (|23I24|) can be used to extract the retarded and advanced component of c. By definition 
c R (t;t') = 5(t-t')c S (t) + 9(t-t')[c > (t;t')-c < (t;t')} 

/>Q<Q 

= a s {t)b s (t')S{t-t') + 6{t-t') / dia R (t;t)[b > (t;t')-b < (i;t')} 

Jo 

poo 

+9{t-t') I di[a>(t;t) - a < (t;t)]b A (i;t'). 



Using the definitions (|21|) and (|22|l to expand the integrals on the r.h.s. of this equation, it is straightforward to show 
that 

c R = a R ■ b R . (25) 

Proceeding along the same lines, one can show that the advanced component is given by c A ~ a A ■ b A . It is worth 
noting that in the expressions for c R and c A no integration along the imaginary track is required. 

Next, we show how to extract the components c) and We first define the Matzubara function fc M (r;r') with 
both the arguments in the interval (0, —i/3): 

fc M (r;r') = k(z = r; z' = t'). 
Let us focus on k) . Without any restrictions we may take t- as the first argument in l|19fl • In this case, we find 

c ] (t;r) = a s (t)b ] (t;T) + dza>(t_; 0)6<(z; r) + / dza<(i_; z)b<(z; r) + / dz a< z)b(z; r). (26) 

Jt + Jo + 

Converting the contour integrals in integrals along the real time axis and along the imaginary track, and taking into 
account the definition in (|21|l 

cl =a R -fcl +a) *b M . (27) 

The relation for can be obtained in a similar way and reads cl" = ■ b A + a M *b^ . Finally, it is straightforward to 
prove that the Matzubara component of c is simply given by c M = a M ★ b M . 

There is another class of identities we want to discuss for completeness. We have seen that the convolution l|18(l of 
two functions belonging to the Keldysh space also belongs to the Keldysh space. The same holds true for the product 

c(z; z') = a(z; z')b(z'] z). 

Omitting the arguments of the functions, one readily finds (for z ^ z') 

c > =a > b < , c < =a < 6 > , J=aW, c^ = a^b\ c M =a M b M . (28) 

The retarded function is then obtained exploiting the identities in l|28|) . We have (for t ^t') 

c R (t; t') = 6(t - t')[a>(t; t')b<(t'; t) - a<(t; t')b>{t'- 1)]. 
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TABLE I: Table of definitions of Keldysh functions and identities for the convolution and the product of two functions in the 
Keldysh space. 



Definition 



c(z; z) = dz a(z; z)b(z; z') 



c[z\ z') — a(z; z')b(z ; z) 



fc > (t;t / ) = k(t+;t'_) 
k<(t;t') = k(t-;f+) 
k R (t;t') = 5{t - t')k s {t) 

+6{t-t')\k > (t;t') -k< (*;*') 
k A (t;t') =S(t-r)k s (t) 

-9(t' — *)[*:>(*;*') - k<(t;t') 
fc 1 (t;r) = fc(t±;r) 
fc r (r,t) = fc(r;t±) 
fc M (r; t') = k(z — r; z' — r') 



c> = a> -b A + a R -b> +aUtf 
c< = a< -b A + a B -b< +a) *tf 



„> 



C = a > fe < 
c< = a < b > 
a R b < +a < b A 
a R b> +a>b A 
a A b< + a<b K 
a A b> +a>b R 
c 1 = Jfe r 
c r = Jb 1 



a ■ o 



c r = a r -& A + a M *6 r 



We may get rid of the 9- function by adding and subtracting a < b < or a > b > to the above relation and rearranging the 
terms. The final result is 

c R = a R b< + a<b A = a R b> + a y b A . 

Similarly one finds c A = a A b < + a < b R = a A b > + a > b R . The time-ordered and anti-time-ordered functions can be 
obtained in a similar way and the reader can look at Table [I] for the complete list of definitions and identities. 

For later purposes, we also consider the case of a Keldysh function k(z; z') multiplied on the left by a scalar 
function l(z). The scalar function is equivalent to the singular part of a function belonging to Keldysh space, 
l(z; z') — l(z)S(z, z'), meaning that f R / A — l M — I and = V = F = 0. Using Tabled one immediately realizes that 
the function / is simply a prefactor: J dz l(z; z)k x (z; z') = l(z)k x (z; z'), where x is one of the Keldysh components 
R, A, ], |~, M). The same is true for f k*(z; z)r(z; z') = k x (z; z')r(z'), where f(z;z') — r(z)S(z, z') and r(z) is a 
scalar function. 



V. THE KADANOFF-BAYM EQUATIONS 

The Green function, as defined in ©, satisfies the equation of motion 

i—G(z;z') = 15(z,z')+h(z)G(z;z')+ [ dzH(z,z)G(z;z'), (29) 
dz J 1 

as well as the adjoint equation 

-i—G(z; z') = 16(z, z') + G(z; z')h(z') + [ dz G(z; z)H(z, z'). (30) 
az' J 1 

The external potential is included in h, while the self-energy X is a functional of the Green function, and describes 
the effects of the electron interaction. The self-energy belongs to Keldysh space and can therefore be written on 
the form S(z,z') = 6(z, z')S 5 (z) + 0(z,z')H > (z ) z') + 0(z' , z)£<(z, z 1 ). The singular part of the self-energy can 
be identified as the Hartree-Fock potential, S 5 (z) = vh(z) + £ x (z). The self-energy obeys the same anti-periodic 
boundary conditions at z = 0_ and z = —i/3 as G. We will discuss self-energy approximations in more detail below. 

Calculating the Green function on the time-contour now consists of two steps: 1) First one has to find the Green 
function for imaginary times, which is equivalent to finding the equilibrium Matzubara Green function G M (r, t'). 
This Green function depends only on the difference between the time-coordinates, and satisfies the KMS boundary 
conditions according to G M (r + i(3,t') — — e l3flN G M (T, r'). Since the self-energy depends on the Green function, 
this amounts to solving the finite-temperature Dyson equation to self-consistency. 2) The Green function with one 
or two time- variables on the real axis can now be found by propagating according to i|29f) and l|3UI) . Starting from 
t = 0, this procedure corresponds to extending the time-contour along the real time-axis. The process is illustrated 
in Fig. |3 Writing out the equations for the components of G using Table [I] we obtain the equations known as the 
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FIG. 2: Propagating the Kadanoff-Baym equations means that one first determines the Green function for time-variables along 
the imaginary track. One then calculates the Green function with one or two variables on an expanding time-contour. 



Kadanoff-Baym equations , 



i-G^(t;t')=h(t)G>(t;t') 



£ R -G§ 



(*;*') 



£^-G A 



(t;t') + 



S 1 *G r 



(31) 



G^-5^ 



(t;f) 



G R -5^ 



(t;t') 



G 1 *£ r 



(32) 



r) = h(t)Gl (t;r) 



and 



G r (r;«)h(t) 



Gl 



sr-G^ 



(t 5 T) + 



(r,t) 



Si *G M 



(*,t), 



(r;t). 



(33) 



(34) 



It is easily seen that if we denote by T the largest of the two time- arguments t and t' , then the right hand sides of 
(|3*T)l and (|32l depend on G^(ti;t 2 ), GI"(ti,£ 2 ) and G"I(£i,t 2 ) for ii,i 2 < T. When propagating the Kadanoff-Baym 
equations one therefore starts at t = t' = 0, with the initial conditions given by G < (0;0) = lim^_>o G M (0; — irf), 
G>(0;0) = lim,^ o G M (-^;0), G^r.O) = G M (r,0) and G1(0,t) = G M (0,r). One then calculates G%i') for 
time-arguments within the expanding square given by t,t' < T. Simultaneously, one calculates G^(t,r) and G^(r, t) 
for t <T. The resulting G then automatically satisfies the KMS boundary conditions. The Kadanoff-Baym equations 
(|31H and l|33|l can both be written in the form 



*-G x (t; z') = h HF (i)G x (t; z') + I x (t; z'), 



(35) 



while (|32|l and (|34|l can be written as the adjoint equations. The term proportional to h HF = h + describes a 
free-particle propagation, while I x is a collision term, which introduces memory effects and dissipation. As can be 
seen from (|31H34|) . the only contribution to I x (0;0) comes from terms containing time-arguments on the imaginary 
axis. These terms therefore contain the effect of initial correlations, since the time-derivative of G would otherwise 
correspond to that of a non-interacting system, i.e., I x (0; 0) = 0. 

An example of a time-propagation is given in Fig. [21 which shows the Green function for an H 2 molecule. In 
this example, the Green function is represented in a basis of Hartree-Fock molecular orbitals, (f\G(z, z')\r l ) = 
4>i(r)Gij( z i z')<f>^{f), where 4>i{r) = (r\4>i)- We have here propagated the Kadanoff-Baym equations using the sec- 
ond Born approximation, illustrated in Fig.0]D. The plots show the imaginary part of the matrix element G^r a (t,t') 
calculated for time- variables within the square t\,t% < T = 20.0 atomic units. In the plot to the left, there is no 
added external potential and the molecule remains in equilibrium. This means that the Green function only depends 
on the difference t 2 —ti, and the oscillations in this time-coordinate has a frequency given by the ionization potential 
of the molecule, in agreement with equilibrium Green function theory |l5j . The figure on the right shows the same 
matrix element, but now in the presence of an additional electric field which is switched on at t = 0. The oscillations 
along the ridge t\ — t 2 can be interpreted as oscillations in the occupation number. 
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FIG. 3: The figures show the Green function ImGj _ (ti,ta), where the matrix indices refer to the groundstate ct 9 Hartree- 
Fock orbital of the molecule. The figure on the left shows the system in equilibrium, while the system on the right has an 
additional electric field, 6(t)Euz. The times t\ and £2 on the axes are given in atomic units. 



VI. CONSERVING APPROXIMATIONS 



In the Dyson-Schwinger equations i|29|) and (|3U[1 . we introduced the electronic self-energy functional S, which 
accounts for the effects of the electron interaction. The self-energy is a functional of the Green function, and will 
have to be approximated in practical calculations. Diagrammatic techniques provide a natural scheme for generating 
approximate self-energies and for systematically improving the approximations. There are no general prescriptions for 
how to select the relevant diagrams, which means that this selection must be guided by physical intuition. There are, 
however, important conservation laws, like the number conservation law or the energy conservation law, that should 
always be obeyed. We will in the following discuss an exact framework for generating such conserving approximations. 

Let us first discuss the conservation laws obeyed by a system of interacting electrons, in an external field given by 
the electrostatic potential v{f,t) and vector potential A(f,t). An important relation between these two quantities is 
provided by the continuity equation 



dt 



n(r,t) + V-j(r,t) = 0. 



(36) 



The density and the current density can be calculated from the Green function using O and (JT3J. Whether these 
quantities will agree with the continuity equation will depend on whether the Green function is obtained from a 
conserving self-energy approximation. If we know the current density we can also calculate the total momentum and 
angular momentum expectation values in the system from the equations 



P{t) = / drj(r,t) 



and 



L(t) = / dff X j(r,t). 



For these two quantities the following relations should be satisfied 



—P(t) = - I dr n(r,t)E{r,t)+j{r,t) x B(P,t) 



dt 



L(t) = - / dr n(f, t)r x E(r, t)+fx (J (f, t) x B(f, t)) 



(37) 

(38) 
(39) 



where E and B are the electric and magnetic fields calculated from 

E(f, t) = Vv(f, t) - d t A(f, t) and B(r, i) = Vx A(f, t). 



(40) 



The equations (|38|l and (|39|l tell us that the change in momentum and angular momentum is equal to the total 
force and total torque on the system. In the absence of external fields these equations express momentum and 
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<n>@ 



(b) 




FIG. 4: Diagrams for the generating functional $, and the corresponding self-energy diagrams. In (a) we have the exchange 
diagram, and (b) the second Born approximation. The diagrams in (c) and (d) belong to the GW approximation and the 
T-matrix approximation respectively. 




angular momentum conservation. Since the right hand sides of l|38H and l)39|l can also directly be calculated from the 
density and the current and therefore from the Green function, we may wonder whether they are satisfied for a given 
approximation to the Green function. 

Finally we will consider the case of energy conservation. Let E(t) = (H{t)) be the energy expectation value of the 
system, then we have 



df 



E(t) 



dfj(f,t) ■ E(r,t) 



(41) 



This equation tells us that the energy change of the system is equal to the work done on the system. The total energy 
is calculated from the Green function using the expression 



G<(t,t')\r) 



(42) 



The question is now whether the energy and the current calculated from an approximate Green function satisfy the 
relation in (|4ip. 

Baym and Kadanoff [l^ . ll4j | showed that conserving approximations follow immediately if the self-energy is obtained 
as the functional derivative, 



S(l;2) 



<5$ 



5G(2;1) 



(43) 



Here, and in the following discussion, we use numbers to denote the contour coordinates, such that 1 = (r*i,^i). A 
functional $ can be constructed, as first shown in a seminal paper by Luttinger and Ward |l2^ . by summing over 
irreducible self-energy diagrams closed with an additional Green function line and multiplied by appropriate numerical 
prefactors, 



*[^=Ei / dl<££W(I;2)G(2; 



1 



(44) 



(k) 

In this summation, £„ denotes a self-energy diagram of n-th order, i.e. containing n interaction lines. The time- 
integrals go along the contour, but the rules for constructing Feynman diagrams are otherwise exactly the same as 
those in the ground-state formalism [Tj| . Notice that the functional derivative in (|43[) may generate other self-energy 
diagrams in addition to those used in the construction of 4> in (|44|l . In Fig. 0] we show some examples of typical i> 
diagrams. Examples of ^-derivable approximations include Hartree-Fock, the second Born approximation, the GW 
approximation and the T-matrix approximation. 

When the Green function is calculated from a conserving approximation, the resulting observables agree with 
the conservation laws of the underlying Hamiltonian, as given in 

GUI), EE), £23), and (S3}. This guarantees the 
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conservation of particles, energy, momentum, and angular momentum. All these conservation laws follow from the 
invariance of $ under specific chan ges in G. We will here only outline the principles of the proofs, without going into 
the details, which can be found in [ll3. IT^j. 1) Number conservation follows from the gauge invariance of $. A gauge 
transformation A(l) — > A(l) + VA(1), where A(r, 0_) = A(r, —ij3) leaves $ unchanged. A consequence of the gauge 
invariance is that a pure gauge cannot induce a change in the density or current. The invariance is therefore closely 
related to the Ward-identities and to the /-sum rule for the density response function 2) Momentum conservation 
follows from the invariance of <I> under spatial translations, r — ► r + R(z). The invariance is a consequence of the 
electron interaction v(l,2) = 5{z\, Z2}/\fi — f%\ being instantaneous and only depending on the difference between 
the spatial coordinates. 3) Angular momentum conservation follows from the invariance of $ under a rotation of the 
spatial coordinates. 4) Energy conservation follows from the invariance of $ when described by an observer using 
a "rubbery clock", measuring time according to the function s(z). The invariance relies on the electron interaction 
being instantaneous. 



VII. NON-INTERACTING ELECTRONS 



In this Section we focus on non-interacting electrons. This is particularly relevant for TDDFT, where the electrons 
are described by the non-interacting Kohn-Sham system. While the Kohn-Sham Green function differs from the true 
Green function, they both produce the same time-dependent density. This is important since the density is not only 
an important observable in, e.g., quantum transport, but also since the density is the central ingredient in TDDFT. 
The use of NEGFs in TDDFT is therefore important due to the relation between u xc and the self-energy. 

For a system of non-interacting electrons H v = and it is straightforward to show that the Green function obeys 
the equations of motion (|29[1 and (|3L)|) . with X = 0. For any z ^ z', the equations of motion can be solved by using 
the evolution operator on the contour, 

S{z,z') =T c {e- i ^' df h ^}, 
which solves i-4-S(z, z') = h(z)S(z,z') and — i-£yS(z, z 1 ) = S(z, z')h(z'). Therefore, any Green function 

G(z; z') = 9(z, z')S(z, 0_)/^S(0_, z') + B{z' , z)S(z, 0_)/<S(0_, z'), (45) 
satisfying the constraint H16[l on the form 

f> - f< = -a, (46) 

is a solution of the I|29I30II . In order to fix the matrix f > or f < we impose the KMS boundary conditions. The 
matrix h(z) = ho for any z on the vertical track, meaning that S(— ij3, 0_) = e _,3ho . Equation (|15fl then implies 
/ < = — e _ ^( h ° _ ^)/ > i and taking into account the constraint l|46|l we conclude that 

f K = e ff(h -M) + l = i/(h0) ' 

where f(u>) = l/[e^^ _M ' + 1] is the Fermi distribution function. The matrix /> takes the form / > = z[/(ho) — 1]. 

The Green function G(z; z') for a system of non-interacting electrons is now completely fixed. Both G > and G < 
depend on the initial distribution function /(ho), as it should according to the discussion of Section ITTT1 Another 
way of writing —iG < is in terms of the eigenstates \<p n ) = |yn(0)) of ho with eigenvalues e n . From the time-evolved 
eigenstate \ip n (t)) = S(£, 0)|^> n ) we can calculate the time-dependent wavefunction ip n (f,t) — {r\ip n (t)). Inserting 
J2 n \fn(Q)){fn{0)\ in the expression for G < we find 

-iG<(r,t-y 7 t') =-i^(r1S(t;Q)^ m )^ m |/ < ^ n )^„|S(Q;i)r') = ^ /(e„)^„(r, t^^, t'), (47) 

rn.n n 

which for t = t' reduces to the time-dependent density matrix. The Green function G > becomes 

-iG>(r, t; ?, = - ^[1 - f{e n )]ip n {?, t'), (48) 

n 

Knowing the greater and lesser Green functions we can also calculate G R,A . By definition we have 

G R (t; t') = 9(t - t') [G> (f; t 1 ) - G< (t; t')] = -i0(t - t')S(t, t'), 
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and similarly 



G A (<; t') = i6(t' - t)S(t, t') = [G R (i' ; i)]t. 



(49) 



In the above expressions the Fermi distribution function has disappeared. The information carried by G R,A is the 
same contained in the one-particle evolution operator. There is no information on how the system is prepared (how 
many particles, how they are distributed, etc). We use this observation to rewrite in terms of G R,A 



G%i') = G R (<;0)G%;0)G A (0;*') 



(50) 



Thus, G^ is completely known once we know how to propagate the one-electron orbitals in time and how they are 
populated before the system is perturbed 0, 0, 0] . We also observe that an analogous relation holds for G^ ■ I" 

Gl(t;i-) =iG R (i;0)Gl(0;r), G*(r;t) = -iG^r; 0)G A (0; r). 

Let us now focus on a special kind of disturbance, namely h(t) = ho + 0(t)hi. In this case 

G R (i; t 1 ) = -i9(t - ^e-^o+hiX*-*') (51) 

depends only on the difference between the time arguments. Let us define the Fourier transform of G R,A from 

G n,A (t; t') = / ^ e -^(*- t ')G RA (^). 

J 2-7T 



with r\ an infinitesimally small positive constant. 



The step function can be written as 9(t — t') = J jjH e J_ i — -. 
Substituting this representation of the ^-function into l|51|l and shifting the to variable one readily finds 

G» = 



h - hi + ir]' 



and therefore G R (ui) is analytic in the upper half plane. On the other hand, from l|49[) it follows that G A (ui) = [G R (u)]^ 
is analytic in the lower half plane. What can we say about the greater and lesser component? Do they also depend only 
on the difference t — t'l The answer to the latter question is negative. Indeed, we recall that they contain information 
on how the system is prepared before hi is switched on. In particular the original eigenstates are eigenstates of h 
and in general are not eigenstates of the Hamiltonian h + hi at positive times. From (|50|l one can see that G>(t; t') 
cannot be expressed only in terms of the time difference t — if . For instance 

G<{t;t') = e- 4(ho+hl ^/(h )e ,; ( ho+hl >*', 

and unless ho and hi commute, it is a function of t and t' separately. 

It is sometimes useful to split h(t) in two parts and treat one of them perturbatively. Let us think, for instance, 
of a system composed of two connected subsystems A + B. In case we know how to calculate the Green function 
of the isolated subsystems A and B, it is convenient to treat the connecting part as a perturbation. Thus, we write 
h(t) = £(t) + V(t), and we define g as the Green function when V = 0. The g is a solution of 



i—- £ (z)} S (z;z') 



16(z,z') 



and of the corresponding adjoint equation of motion. Furthermore, the Green function g obeys the KMS boundary 
conditions. With these we can use g to convert the equations of motion for G into integral equations 



G(z;z')=g(z;z')+ f dz g(z; z)V(z)G(z; z') = g(z; z') + f dz G(z; z)V(z)g(z; z'); 



(52) 



the integral on z is along the generalized Keldysh contour of Fig. ^ One can easily check that this G satisfies both 
(|29|l and G also obeys the KMS boundary conditions since the integral equation is defined on the contour of 

Fig. [[] 

In order to get some familiarity with the above perturbation scheme, we consider explicitly the system A+B already 
mentioned. We partition the one-electron Hilbert space in states of the subsystem A and states of the subsystem B. 
The "unperturbed" system is described by £, while the connecting part by V and 



s = 


£aa 


v = 


Vab 




S BB 




Vba 
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Taking into account that g has no off-diagonal matrix elements, the Green function projected on one of the two 
subsystems, e.g., G BB , is 

G BB (z;z')=g BB (z;z') + / dz g BB (z; z)Vba(z)G A b{z; z') 



and 



G A b{z\z') = I dz %aa{z;z)Vab{z)Gbb{z\z'). 
Substituting this latter equation into the first one, we obtain a closed equation for G BB : 

G BB (z;z')=g BB (z;z') + dzdz'g BB (z; z)T, BB (z; z')G BB (z'; z'), (53) 

J 7 

with 

S b _b(z; z ) = Vba(z)eaa(z; z')Vab(z) 

the embedding self-energy. The retarded and advanced component can now be easily computed. With the help of 
Tabled one finds 

piR,A _ R,A , R.A V R.A p,R,A 
"BB — &BB &BB ' BB ' BB ' 



Next, we have to compute the lesser or greater component. As for the retarded and advanced components, this 
can be done starting from (|53|l . The reader can soon realize that the calculation is rather complicated, due to the 
mixing of pure real-time functions with function having one real time argument and one imaginary time argument, 
see Tabled Below, we use (|50f) as a feasible short-cut. A closed equation for the retarded and advanced component 
has been already obtained. Thus, we simply need an equation for G^(0;0). Let us focus on the lesser component 
G < (0;0) = if < . Assuming that the Hamiltonian ho is hermitian, the matrix (ui — ho) -1 has poles at frequencies 
equal to the eigenvalues of ho. These poles are all on the real frequency axis, and we can therefore write 

G<(0;0)=i/(h ) = / g/(C)-i-, (54) 

where the contour 7 encloses the real frequency axis. 



VIII. ACTION FUNCTIONAL AND TDDFT 



We define the action functional 

A = ilnTr |e /3 ^S'(-i/3; 0)} , (55) 

where the evolution operator S is the same as defined in . The action functional is a tool for generating equations 
of motion, and is not interesting per se. Nevertheless, one should notice that the action, as defined in (|55|l has a 
numerical value equal to ihiZ, where Z is the thermodynamic partition function. 

It is easy to show that if we make a perturbation 5V(z) in the Hamiltonian, the change in the evolution operator 
is given by 

i—5S(z; z') = 5V(z)S(z; z') + H(z)5S(z; z'). (56) 
dz 

A similar equation for the dependence on z 1 , and the boundary condition SS(z; z) — gives 

5S(z;z') = -i [ dzS{z-z)8V{z)S{z,z'). (57) 
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We stress that the time-coordinates are on a contour going from to — i(3. The variation in, e.g., V(t+) is therefore 
independent of the variation in V(t-). If we let 5V(z) = J dr5v(r,z)h(r), a combination of (|55|l and (|57|l yields 
[compare to (@J] the expectation values of the density, 

5A _ I 5 xr [ C 0^S( j<3-0)} 

5v(r,z) Xr |e^5(-i/9;0)} 6v (^ z ) 1 J 

Tr f e^S(-i/3; 0)5(0; z)n(f)S(z, 0) } 

1 — i -=n{r,z). (58) 



Tr 



l^/^^-ijS; 0)} 



A physical potential is the same on the positive and on the negative branch of the contour, and the same is true for 
the corresponding time-dependent density, n(r, t) = n(r, t±). A density response function defined for time-arguments 
on the contour is found by taking the functional derivative of the density with respect to the external potential. Using 
the compact notation 1 = (ri,zi), the response function is written 

^ = W> = W)W) =xi2 '- iy (59) 

This response function is symmetric in the space and time-contour coordinates. We again stress that the variations 
in the potentials at t + and t_ are independent. If, however, one uses this response function to calculate the density 
response to an actual physical perturbing electric field, we obtain 

5n(r, t) = Sn(r, t±) = J dz J drx(r,t±;rz')5v(r,z'), (60) 

where 7 indicates an integral along the contour. In this expression, the perturbing potential (as well as the induced 
density response) is independent of whether it is located on the positive or negative branch, i.e. Sv(r* , = Sv(r' , t'). 
We consider a perturbation of a system initially in equilibrium, which means that 5v{r l \t') =/= only for t' > 0, and 
we can therefore ignore the integral along the imaginary track of the time-contour. The contour integral then consists 
of two parts: 1) First an integral from t' = Q to t' = t, in which \ — X > : an d 2) an integral from t' = t to t' = 0, 
where x = X < - Writing out the contour integral in H60[) explicitly then gives 

Sn(r,t)= [ dt' [ df'[x > (rt;r / t')-x < (rt;rt')]5v(f',t')= [ dt' f dr 1 X R (rt;r't')5v(f' ,t'). (61) 



The response to a perturbing field is therefore given by the retarded response function, while x(l> 2) defined on the 
contour is symmetric in (1 <-> 2). 

If we now consider a system of non-interacting electrons in some external potential v Sl we can similarly define 
a non- interacting action- functional A s . The steps above can be repeated to calculate the non- interacting response 
function. The derivation is straightforward, and gives 

Xs(l;2) = ^|^ = -iG 8 (l;2 )Gs( 2;l ) . (62) 

The non- interacting Green function G s has the form given in l|45|) . 147|) and (|48|l . The retarded response-function is 

Xf(**i,*i;»*2,*2) - -i9(ti-t 2 ) [G> {r 1 M\r 2 M)G<{r 2 ^r 1 M) - Gf{r 1 M-,r 2 M)G> {r 2 M]riM)\ 

= i^2[f{£m) - f(e n )]<Pn(r 1 ,t 1 )ip* m (r 1 ,t 1 )tp m (r2,t2)tPn(r2,t 2 ), (63) 

where we have used l|47|l and (|48|l in the last step. 

Having defined the action functional for both the interacting and the non-interacting systems, we now make a 
Legendre transform, and define 



A[n] = -A[v] + J d(l)n(l)u(l), (64) 
which has the property that 8A[n]/Sn(l) = v(l). Similarly, we define the action functional 

A s [n} = -A s [v s }+ [ d(l)n(l)v s (l). (65) 
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with the property SA s [n]/8n(l) — v s (l). The Legendre transforms assume the existence of a one-to-one correspondence 
between the density and the potential. From these action functionals, we now define the exchange-correlation part to 
be 

A xc [n] = A s [n] - A[n] - \ [ d(12)5(z 1 ,z 2 ) 1 &^\ . (66) 

Taking the functional derivative with respect to the density gives 

v s [n]{l) = w(l) + uh(1) + v xc [n](l) (67) 

where ujj(l) is the Hartree potential and v xc (l) — 8A xc /8n(X). Again, for time-arguments on the real axis, these 
potentials are independent of whether the time is on the positive or the negative branch. If we, however, want to 
calculate the response function from the action functional, then it is indeed important which part of the contour the 
time-arguments are located on. 

We already described how to define response function on the contour, both in the interacting (|59|l and the non- 
interacting (|62[1 case. Given the exact Kohn-Sham potential, the TDDFT response function should give exactly the 
same density change as the exact response function, 



Sn(l) = / d(2) X (l; 2)8v(2) = / d(2) Xs (l; 2)Sv a (2). (68) 



The change in the Kohn-Sham potential is given by 



= 6v(l) + J d(2) /hxc(1; 2)5n(2) - 8v(l) + J d(23) / Hxc (l; 2) X (2; 3)5v(S) (69) 
where /hxc(1; 2) = 8(zi, z 2 )/\r\ — r 2 \ + 8v xc (l) / 8n(2) . Inserted in l|68|) . we obtain 

X (l; 2) = Xs (l; 2) + / d(34) Xs (l; 3)/ Hx c(3; 4) X (4; 2). (70) 



This is the response function defined for time-arguments on the contour. If we want to calculate the response induced 
by a perturbing potential, the density change will be given by the retarded response function. Using Table [IJ we can 
just write down 

X R (ri,ti;r 2 ,t 2 ) = xf (n, h; fb, t 2 ) + / dtzdtidrzdr^ xf {f\, t x ; r 3 , t 3 )f§ xc (r 3 , t 3 ; r 4 , t 4 )x R (r±, U; r 2 , t 2 ). (71) 



The time-integrals in the last expression go from to oo. As expected, only the retarded functions are involved in 
this expression. We stress the important result that while the function /h X c(1,2) is symmetric under the coordinate- 
permutation (1 «-> 2), it is the retarded function 

fi KC (ruh;r 2 ,t 2 ) = ^iiM + fK(? u h; r 2 ,t 2 ) (72) 
\n - r 2 \ 

which is used when calculating the response to a perturbing potential. 

IX. EXAMPLE: TIME-DEPENDENT OEP 

We will close this section by discussing the time-dependent optimized effective potential (TDOEP) method in the 
exchange-only approximation. This is a useful example of how to use functions on the Keldysh contour. While the 
TDOEP equations can be derived from an action functional, we will here use the time-dependent Sham-Schluter 
equations as starting point [17j. This equation is derived by employing a Kohn-Sham Green function, G s (l, 2) which 
satisfies the equation of motion 

i ~Iz~ + ~2 ~ v *(^> z ^) \ Gsi^i, zi; r 2 , z 2 ) = 8(z 1 , z 2 )S(fi - f 2 ), (73) 
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as well as the adjoint equation. The Kohn-Sham Green function is given by (|47|l and (|48|l in terms of the time- 
dependent Kohn-Sham orbitals. Comparing (|73[) to the Dyson-Schwinger equation (|29[) . we see that we can write an 
integral equation for the interacting Green function in terms of the Kohn-Sham quantities, 

G(l; 2) = G s (l; 2) + J d(12) G s (l; 1) [S(l; 2) + 5(1, 2>(1) - v s (l)}} G(2; 2). (74) 

It is important to keep in mind that this integral equation for G(l; 2) differs from the differential equations l|29|) and 
(1301) in the sense that we have imposed the boundary conditions of G s on G in l|74|l . This means that if G s (l;2) 
satisfies the KMS boundary conditions l)15[l. then so will G(l;2). 

If we now assume that for any density n(l) = —iG(l; 1 + ) there is a potential v s (l) such that n(l) = — iG s (l; 1 + ), 
we obtain the time-dependent Sham-Schluter equation, 



d(12) G 8 (l; 1)S(1; 2)G(2; 1) = / d(I)G.(l; I) [«„(!) - u(l)]G(l; 1). (75) 



This equation is formally correct, but not useful in practice since solving it would involve first calculating the nonequi- 
librium Green function. Instead, one sets G — G s and E[G] = S[G S ]. For a given self-energy functional, we then 
have an integral equation for the Kohn-Sham equation. This equation is known as the time-dependent OEP equation. 
Defining £ = vh + E xc and v s = v + «h + w xc , the TDOEP equation can be written 

/'d(12)G s (l;l)E xc [G s ](l;2)G 8 (2;l) = J d(T)G s (l; 1> XC (T)G S (1; 1). (76) 

In the simplest approximation, S xc is given by the exchange-only self-energy of Fig. 

Ex(l; 2) = iG< (1; 2)t>(l, 2) = - £ n,^-(l)0*(2) W (l, 2) (77) 

j 

where rij is the occupation number. This approximation leads to what is known as the exchange-only TDOEP 
equations [iH. IT9I l2(i | . Since the exchange self-energy S x is local in time, there is only one time-integration in H76fl . 
The x-only solution for the potential will be denoted v x . With the notation S(3; 4) = E x (r 3 i 3 ; f^t^) — 8(fs — ri)v x (r3t3) 
we obtain from l|76|l 

11 dt 3 J dr 3 dr A [G<(1;3)S(3;4)G>(4;1) - G>(1;3)E(3;4)G<(4; 1)" 

dt 3 / dr 3 ^4G](l;3)E(3;4)G[(4;l). (78) 



Let us first work out the last term which describes a time-integral from to —i[3. On this part of the contour, the 
Kohn-Sham Hamiltonian is time- independent, with v x (f, 0) = v x (r), and (fi(r,t) = <pi(r) cxp (— idt). Since E x is 
time-independent on this part of the contour, we can integrate 

-if! e P(ei-ek) _ 1 

dt 3 G](l;r 3 ,t 3 )Gl(r 4 ,t 3 ;l) = -i Vn,(l - n fc )^(l)^(r 3 )^(r 4 )^(l) (79) 

— Si ~ £k 

i,k 

If we then use 7ij(l — nfc)(e /3 ' £i ~ £fe ' ) — 1) = nk — ni and define the function u x j by 



^•(i) = -^E^/ 



d2^(2)^(2K(lV(l,2) (80) 



we obtain from fT^jl and fTTjl 

„ f wi,, ^„r, A ^ f v-^ ^(^2)^(^2) 



/ 4 di 3 / dr 3 dr 4 G](l,3)E(3,4)G[(4,l) = - / rff 2 ^ Wj ^ ^^ ^(1)^(1) K 3 (r 2 ) - ^(r 2 )] 



c.c. 



(81) 
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The integral along the real axis on the lhs of (|78|l can similarly be evaluated. Collecting our results we obtain the 
OEP equations on the same form as in [20l| . 



= i 

3 k^j 



r dt 2 f dr 2 [v K (2)-Ux, j (2)}p j (l)<p*(2)<p* k (l)<p k (2)+c.c. 

3 k*j J ° J 

.££ n .w(iM(i) 



£j — £k 

3 fc/j J 



dr 2 <p* (f 2 ) [v x (f 2 ) - u x j (f 2 )] tp k (r 2 ). (82) 



The last term represents the initial conditions, expressing that the system is in thermal equilibrium at t = 0. The 
equations have exactly the same form if the initial condition is specified at some other initial time to ■ ln the case that 
we let t n — * -co, the term due to the initial conditions vanish and the remaining expression equals the one given in 
\17l llSl H9| . The OEP-equations l|82|l in the so-called KLI-approximation have been successfully used by Ullrich et 
al.|19l| to calculate properties of atoms in strong laser fields. 
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